library(foreign)
abgss = read.dta("~/Documents/abgss2.dta")

install.packages("glmx", repos="http://R-Forge.R-project.org")
install.packages("lmtest")
library(lmtest)
library(glmx)
attach(abgss)

##I found a version of heterskedastic probit analysis code for R that was prepared for the original article, but I couldn't get the code to work. 
#Here are some of my tinkerings with the code. 

ab1 =glm(abany ~ reliten + attend + male + prot + cath + jew + era,
         family="binomial"(link=probit))
sigma = cbind(abproct, abconct, abboth, abinfo)
ivs = cbind(reliten, attend, male, prot, cath, jew, era)
summary(ab1)

hetprob(ivs, sigma, abany) #This will be attached as a separate program.

#### glmx package

health <- hetglm(abhlth ~ black + male + cath + reliten + attend + erameans + era  |
  abproct*abconct  + abinfo + abimp+ abfirm, data = abgss)
summary(health)
#attempt similar version of the model using logit, no assumption of heteroskedasticity

health.log = glm(abhlth ~ black + male + cath + reliten + attend + erameans + era  +
  abproct*abconct  + abinfo + abimp+ abfirm, data = abgss, family=binomial(link = "logit"))
summary(health.log)

#attempt similar version of the model using probit, ignores heteroskedasticity

health.pro = glm(abhlth ~ black + male + cath + reliten + attend + erameans + era  +
  abproct*abconct  + abinfo + abimp+ abfirm, data = abgss, family=binomial(link = "probit"))
summary(health.pro)

library(apsrtable)
summary(health)


poor <- hetglm(abpoor ~ black + male + cath + reliten + attend + erameans + era  |
  abproct*abconct  + abinfo + abimp+ abfirm, data = abgss)
summary(poor)
rape = hetglm(abrape ~ black + male + cath + reliten + attend + erameans + era  |
  abproct*abconct  + abinfo + abimp+ abfirm, data = abgss)
summary(rape)
single = hetglm(absingle ~ black + male + cath + reliten + attend + erameans + era  |
  abproct*abconct  + abinfo + abimp+ abfirm, data = abgss)
summary(single)
nomore = hetglm(abnomore ~ black + male + cath + reliten + attend + erameans + era  |
  abproct*abconct  + abinfo + abimp+ abfirm, data = abgss)
summary(nomore)
defect = hetglm(abdefect ~ black + male + cath + reliten + attend + erameans + era  |
  abproct*abconct  + abinfo + abimp+ abfirm, data = abgss)
summary(defect)
any = hetglm(abany ~ black + male + cath + reliten + attend + erameans + era  |
  abproct*abconct  + abinfo + abimp+ abfirm, data = abgss)
summary(any)


table <- sapply(names(abgss)[1:7], function(x) {
  f <- as.formula(paste(x,
                        "~ black + male + cath + reliten + attend + erameans + era",
                        "| abproct*abconct  + abinfo + abimp+ abfirm"))
  f0 <- as.formula(paste(x, "~ 1"))
  m <- hetglm(f, data = abgss)
  m0 <- hetglm(f0, data = model.frame(m))
  c(Percent_yes = as.vector(100 * prop.table(table(abgss[[x]]))["1"]),
    coef(m)[c(1:10, 14, 11:13)],
    Heteroskedasticity = as.vector(summary(m)$lrtest[1]),
    N = nobs(m),
    Goodness_of_fit = 2 * as.vector(logLik(m) - logLik(m0))
  )
})
round(table, digits = 2)
xtable(tab1)

objects(health)


#Estimating True Values
provec = seq(min(abgss$abproct), max(abgss$abproct), length.out=100)

predvec.health <- NULL
for(i in 1:length(provec)){
  values <- c(1, 0, 1, 0, mean(abgss$reliten, na.rm=T), mean(abgss$attend, na.rm=T), 1, median(abgss$era, na.rm=T), 
              provec[i], median(abgss$abconct, na.rm=T), median(abgss$abinfo, na.rm=T), median(abgss$abimp, na.rm=T),
              median(abgss$abfirm, na.rm=T), median(abgss$abproct*abgss$abconct, na.rm=T))
  pred.health <- coef.hetglm(health) %*% values
  predvec.health <- c(predvec.health, pred.health)
}

predprob.health <- 1/(1+exp(-predvec.health))

predvec.poor <- NULL
for(i in 1:length(provec)){
  values <- c(1, 0, 1, 0, mean(abgss$reliten, na.rm=T), mean(abgss$attend, na.rm=T), 1, median(abgss$era, na.rm=T), 
              provec[i], median(abgss$abconct, na.rm=T), median(abgss$abinfo, na.rm=T), median(abgss$abimp, na.rm=T),
              median(abgss$abfirm, na.rm=T), median(abgss$abproct*abgss$abconct, na.rm=T))
  pred.poor <- coef.hetglm(poor) %*% values
  predvec.poor <- c(predvec.poor, pred.poor)
}

predprob.poor <- 1/(1+exp(-predvec.poor))

predvec.rape <- NULL
for(i in 1:length(provec)){
  values <- c(1, 0, 1, 0, mean(abgss$reliten, na.rm=T), mean(abgss$attend, na.rm=T), 1, median(abgss$era, na.rm=T), 
              provec[i], median(abgss$abconct, na.rm=T), median(abgss$abinfo, na.rm=T), median(abgss$abimp, na.rm=T),
              median(abgss$abfirm, na.rm=T), median(abgss$abproct*abgss$abconct, na.rm=T))
  pred.rape <- coef.hetglm(rape) %*% values
  predvec.rape <- c(predvec.rape, pred.rape)
}

predprob.rape <- 1/(1+exp(-predvec.rape))


predvec.single <- NULL
for(i in 1:length(provec)){
  values <- c(1, 0, 1, 0, mean(abgss$reliten, na.rm=T), mean(abgss$attend, na.rm=T), 1, median(abgss$era, na.rm=T), 
              provec[i], median(abgss$abconct, na.rm=T), median(abgss$abinfo, na.rm=T), median(abgss$abimp, na.rm=T),
              median(abgss$abfirm, na.rm=T), median(abgss$abproct*abgss$abconct, na.rm=T))
  pred.single <- coef.hetglm(single) %*% values
  predvec.single <- c(predvec.single, pred.single)
}

predprob.single <- 1/(1+exp(-predvec.single))


predvec.nomore <- NULL
for(i in 1:length(provec)){
  values <- c(1, 0, 1, 0, mean(abgss$reliten, na.rm=T), mean(abgss$attend, na.rm=T), 1, median(abgss$era, na.rm=T), 
              provec[i], median(abgss$abconct, na.rm=T), median(abgss$abinfo, na.rm=T), median(abgss$abimp, na.rm=T),
              median(abgss$abfirm, na.rm=T), median(abgss$abproct*abgss$abconct, na.rm=T))
  pred.nomore <- coef.hetglm(nomore) %*% values
  predvec.nomore <- c(predvec.nomore, pred.nomore)
}

predprob.nomore <- 1/(1+exp(-predvec.nomore))

predvec.defect <- NULL
for(i in 1:length(provec)){
  values <- c(1, 0, 1, 0, mean(abgss$reliten, na.rm=T), mean(abgss$attend, na.rm=T), 1, median(abgss$era, na.rm=T), 
              provec[i], median(abgss$abconct, na.rm=T), median(abgss$abinfo, na.rm=T), median(abgss$abimp, na.rm=T),
              median(abgss$abfirm, na.rm=T), median(abgss$abproct*abgss$abconct, na.rm=T))
  pred.defect <- coef.hetglm(defect) %*% values
  predvec.defect <- c(predvec.defect, pred.defect)
}

predprob.defect <- 1/(1+exp(-predvec.defect))


predvec.any <- NULL
for(i in 1:length(provec)){
  values <- c(1, 0, 1, 0, mean(abgss$reliten, na.rm=T), mean(abgss$attend, na.rm=T), 1, median(abgss$era, na.rm=T), 
              provec[i], median(abgss$abconct, na.rm=T), median(abgss$abinfo, na.rm=T), median(abgss$abimp, na.rm=T),
              median(abgss$abfirm, na.rm=T), median(abgss$abproct*abgss$abconct, na.rm=T))
  pred.any <- coef.hetglm(any) %*% values
  predvec.any <- c(predvec.any, pred.any)
}

predprob.any <- 1/(1+exp(-predvec.any))

plot(provec, predprob.health, type="l", ylim=c(0,1), xlab="Pro Abortion Arguments", ylab="Probability of Support", main="Probability of Supporting Different Abortion Scenarios")
points(provec, predprob.poor, type="l", col="red")
points(provec, predprob.rape, type="l", col="blue")
points(provec, predprob.single, type="l", col="slategray")
points(provec, predprob.nomore, type="l", col = "green")
points(provec, predprob.defect, type="l", col= "midnightblue")
points(provec, predprob.any, type="l", col="skyblue")
legend("bottomleft",c("Health","Poor", "Rape", "Single", "No More", "Defect", "Any"), 
       col=c("black","red", "blue", "slategray", "green", "midnightblue", "skyblue"), lty=c(1,1), cex=.4)


